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We study the n = 2 Renyi entanglement entropy of the triangular quantum dimer model via Monte 
Carlo sampling of Rokhsar-Kivelson(RK)-like ground state wavefunctions. Using the construction 
proposed by Kitaev and Preskill [Phys. Rev. Lett. 96, 110404 (2006)] and an adaptation of 
the Monte Carlo algorithm described in [Phys. Rev. Lett. 104, 157201 (2010)], we compute 
the topological entanglement entropy (TEE) at the RK point 7 = (1.001 =b .003) In 2 confirming 
earlier results. Additionally, we compute the TEE of the ground state of a generalized RK-like 
Hamiltonian and demonstrate the universality of TEE over a wide range of parameter values within 
a topologically ordered phase approaching a quantum phase transition. For systems sizes that are 
accessible numerically, we find that the quantization of TEE depends sensitively on correlations. 
We characterize corner contributions to the entanglement entropy and show that these are well 
described by shifts proportional to the number and types of corners in the bipartition. 



I. INTRODUCTION 

Quantum liquid phases of matter that do not break 
conventional symmetries can have "hidden" non-local 
quantum orders. Such quantum liquids are ordered quan- 
tum phases that are not described by a local order param- 
eter. Topologically ordered phases in particular, are of 
great interest because of their potential to form the basis 
of a physically fault-tolerant quantum computer. 1 There 
is therefore a strong incentive to realize such phases in 
experimental systems as well as to identify theoretical 
models which posses topologically ordered phases. 

However, the lack of a local order parameter inhibits 
the identification of topological ordered phases in theo- 
retical models. Kitaev & Preskill 2 and Levin & Wen 3 
identified a sub-leading negative constant term in the bi- 
partite entanglement entropy, the topological entangle- 
ment entropy (TEE), which allows for the identification 
and classification of topological order. Constant sublead- 
ing terms can arise in other contexts including critical 
systems from Goldstone modes in symmetry broken 
states and from corners in non-smooth bipartition as 
seen in integer quantum hall wavefunctions. 9 

Lattice models with hard local constraints, such as 
quantum dimer and loop models, possess quantum liq- 
uid ground state phases, including topological phases. In 
particular, the hard-core quantum dimer model on the 
triangular lattice (TQDM) has a Z2 topologically ordered 
dimer liquid phasef^EU Since topological phases gener- 
ally arise in strongly interacting systems which are not 
always tractable by analytic methods, numerical studies 
of these models are often necessary. Lanczos diagonaliza- 
tion may be used to compute the bipartite-entanglement 
entropy in small systems.^ However, computing the sub- 
leading term in the entanglement requires using moder- 
ately large systems which are not generally accessible via 
Lanczos diagonalization. 

At the RK point the TEE of the TQDM has been 



computed using Pfafnan (Kasteleyn) methods with high 
(10 -9 ) numerical accuracy,* 12 * 13 * and away from the RK 
point using exact diagonalization on small lattices 
Recent work by Hastings et al™ has demonstrated a 
method for computing the Renyi entanglement entropy 
via Monte Carlo methods. This is attractive, since these 
techniques generally allow for the study of moderately 
large systems. 

Demonstrating the universality of the TEE within a 
topological phase and its behavior across phase tran- 
sitions, is an area of ongoing research. Temperature 
induced transitions have been explored in the work of 
Isakov et a/P^ The behavior of TEE approaching a quan- 
tum phase transition was previously studied by Stephan 
et al by interpolating between the triangular to the 
square lattice dimer models P 

Here we adapt the method of Hastings el aP^ to 
the TQDM . We confirm the previous results for the 
TEE™ and characterize constant contributions due to 
corner effects at the RK point which may compete with 
the TEE. Additionally we compute the TEE of a "gen- 
eralized" RK wave-function, using the model of Trous- 
selet et al ™ and show the evolution of TEE approach- 
ing a first order quantum phase transition. The results 
strongly suggest the universal nature of the TEE inside 
the dimer liquid phase in the thermodynamic limit, al- 
though correlations are found to limit convergence in fi- 
nite systems. 



A. Triangular Quantum Dimer Model 

The fully-packed hard-core dimer model is defined on 
a lattice with degrees of freedom labeled by the occupa- 
tion of dimers on links, and the constraint that exactly 
one dimer must touch each vertex. The Hilbert space is 
comprised of the fully packed dimer coverings on the lat- 
tice satisfying the vertex constraint (= |{C})). Different 
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dimer coverings are defined to be orthogonal. Rokhsar 
and Kivelson first introduced this model on the square 
latticed It was subsequently generalized by Moessner 
and Sondhi^to the triangular lattice, where indications 
of a Z 2 topologically ordered ground state emerged The 
Hamiltonian for the TQDM is: 

H = J2 + h.c)+v(\&)(&\ + |£7)<£7|) (1) 

P 

where p labels all minimal Rhombus plaquettes on the 
triangular lattice, and the kinetic (t) term which flips 
parallel dimers around a plaquette is the minimal dimer 
hopping that respects the vertex constraint (here t > 
always). The v term acts as a potential energy between 
parallel dimers. At the RK point, defined by v = t, the 
ground state can be written as 

c vz 

where the sum is over configurations reachable by plaque- 
tte flips (Z is the number of elements in the sum). On 
the torus, plaquette flips conserve two parities that are 
defined by counting the occupation of dimers intersect- 
ing the two non-contractible loops of the torus. Dimer 
coverings are split into four topological sectors defined by 
these parities, such that local rearrangements of dimers 
cannot connect configurations in two different topologi- 
cal sectors. Plaquette flips on the triangular lattice are 
believed to be nearly ergodic within a topological sector, 
with the exception of 12 symmetry related "staggered" 
configurations that have no flippable plaquettes. There- 
fore four distinct ground states are defined by these par- 
ities which we label ft = (0,0), (1,0), (0, 1), (1, 1) (0 for 
even parity). This topological degeneracy is a character- 
istic of the topological order present in the ground state 
at the RK point. 

On the triangular lattice the topologically ordered 
dimer liquid phase persists for a finite reg ion b elow the 
RK point, in the range 0.86 ~ v/t < lpEHl Outside 
of this region (v/t < 0.86 and v/t > 1), the ground 
state is one of several symmetry broken ordered crys- 
talline phases.^ 

The RK wave function can be generalized to weighted 
superpositions of dimer configurations {C}, 

i*> = £^ e - E(c) |c>, ( 3 ) 

c ^ Z 

with Z = J2e~ 2E(C) , and E(C) is the "classical" enerev 

c 

of C. Such a generalized RK wave function is the ex- 
act zero energy ground state of a corresponding RK-like 
local Hamiltonian. 19 In Sec. |III^ >, we compute the TEE 
of a generalized RK wave function which was previously 
studied in Refs.H9land l20l where it is seen to interpolate 
between the topologically ordered phase and a symmetry 
broken phase. 



B. Entanglement Entropy and Topological Phases 

Bipartite entanglement entropy has emerged as a pow- 
erful probe of quantum systems. The bipartite entangle- 
ment entropy of a pure state is defined by specifying 
a bipartition of the lattice into a region A and its com- 
pliment B. The von Neumann entropy is defined as 

S(p A ) = -Tip A \np A (4) 
and the Renyi entropy is defined as 

S n { PA ) = j^—Trp n A (5) 
1 — n 

where p A = Trg|\I/)(\I/| is the reduced density matrix 
of A. The Renyi entropy reduces to the von Neumann 
entropy in the limit n —> 1 and both are symmetric under 
exchange of A and B, S n (p A ) — S n (pB)- Ground states 
of local Hamiltonians are known to exhibit a boundary 
law scaling in region size (i.e., in 2D the scaling is with 
the perimeter length) although critical fermions are a 
notable exception to this ruleP^ In two dimensions this 
scaling can generically be written as 

S( PA ) = aL A +l3 log(L A /a) + C + 0(1/ L A ) (6) 

where the leading term is proportional to the perimeter 
L A , and a is a non- universal constant. The logarithmic 
term appears in certain quantum critical theories, but for 
gapped phases, it is expected that f3 = 0. The constant 
term Co has been shown to arise in critical phases as well 
as topologically ordered phases PtZl 

For topological phases, there is a universal, negative, 
constant subleading term, the topological entanglement 
entropy: —7 G Co (with 7, also referred to as "ytopo, 
> 0) Topological phases may be des cribe d by an effec- 
tive topological quantum field theoryf^I] such theories 
are categorized by the so-called total quantum dimen- 
sion, D. For conventional ordered phases D = 1 and for 
topologically ordered phases D > 1. The TEE is given 
by: 

7 = lnL> (7) 

and therefore is an indication of topological order (7 = 
for conventional phases). 

Physically the origin of this term can be seen by con- 
sidering string-net wavefunctions as an effective theory 
of topological order!^ The non-local order encoded in a 
topologically ordered phase can be understood in terms 
of effective loop or string-net degrees of freedom describ- 
ing the wavefunction. Specifically for discrete gauge the- 
ories, wavefunctions are comprised of different types of 
non-branching loops with (counting the absence of a loop 
as one type), the relation: the number of types of loops 
= to the elements of the group = the total quantum di- 
mension, D. Then as a direct consequence of the fact 
that each type of loop must enter and exit the boundary 
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an even number of times, the effective degrees of free- 
dom crossing the bipartition boundary as probed by the 
entanglement entropy, is corrected by a factor of 1/D. 
This is responsible for the reduction of the entanglement 
entropy scaling by lnD. 

The dimer model belongs to the Z 2 class^l Q f topo- 
logically ordered phases 2 ^ 2 ^ (the effect ive loop degrees 
of freedom are so-called transition loops ] 20 * 25 * 26 ! see sup- 
plementary materiaP^). Therefore for the dimer model 
and other Z 2 topologically ordered phases 7 = In 2. Fur- 
thermore it has been shownj 2 ^ that 7 is independent of 
the Renyi parameter n, so that any Renyi entanglement 
entropy can be used to compute the quantity 7. 



In Sec. IV we show that for non-smooth bipartitions 
on a lattice, there can be non-universal constant con- 
tributions to Co. We will split Co into universal and 
non-universal parts by writing Co = —7 + ft. Our re- 
sults are consistent with a non-universal term k of the 
form n = where is the number of corners of 

type i. As discussed in Sec. |IV[ the origin of these corner 
terms can be thought of coming from a substitution in a 
generalized linear scaling Sa ~ Y^i where boundary 
vertices i contribute different constants <%•. 



II. MONTE CARLO SAMPLING FOR 
ENTANGLEMENT ENTROPY 

In Ref. 14 Hastings et al. describe a SWAP algorithm 
to compute S 2 (pa) y i a Monte Carlo simulations (in the 
following, Sa will always be taken to mean S 2 (pa))- In 
the current work we will be considering generalized RK 
points, characterized by wavefunctions that are explicitly 
written as a weighting of configurations, therefore we are 
able to compute expectation values of estimators using 
classical Monte Carlo sampling of the wavefunction. 

To estimate the entanglement entropy, following 
Ref. [14j we define a new "doubled" system as two non- 
interacting independent copies of the original, labeled 1 
and 2. Each corresponding copy has the identical bipar- 
tition A and B so that the Hilbert space of the doubled 
system is a tensor product of the two copies with a state 
labeled by degrees of freedom in Ai, B\ and A 2l B 2 re- 
spectively. Then Sa is related to the expectation of the 
SWAP^ operator defined on the doubled Hilbert space 
by its action in swapping the degrees of freedom in A: 

SWAPaIAx^x) ® \A 2 B 2 ) = \A 2 Bi) <g> |AiS 2 ). (8) 

Hastings et al. showed that Trp\ = (SWAP^), and 
therefore S A = - ln(SWAP A ). 

Taking C to represent a "doubled" dimer covering, the 
matrix elements of the SWAP operator are 



(C'\SWAP a \C) = 5c>cJ(C\a), 



(9) 



where Ca is the configuration resulting from swapping C 
over region A, and S(C\a) is one if the swapped config- 
urations do not violate the hard-core dimer constraint, 




FIG. 1. (Color Online) Regions A' and A as described in the 
text. Note A' also includes the blue links. Shaded plaquettes 
are examples of constrained plaquettes for the case where A 
is kept flippable. 



and zero otherwise. We can write the expectation value 
of SWAP as weighted sum over configurations C: 

(SWAP A ) = \ e- E ^- E ^ c \C\SWAP A \C) 



c,c 



4£. 



-E(C A )-E 



{C) S(C\ a ) 



(10) 



= ^e-™^|,)n(c), 



where Z = Y,c exp(-2E(C)), 11(C) = exp{-2E(C))/Z 
can be viewed as a probability distribution, and 
AE(C\ A ) = E{C A ) ~ E(C). We see then that by clas- 
sical Monte Carlo sampling of 11(C), we can compute 
the expectation value of SWAP by using the estimator 
e~ AE ^S(C\ A )- For the RK ground state, the expecta- 
tion value of the swap operator is simply the fraction of 
the dimer configurations that are A-swappable. 

As a direct consequence of the perimeter law scaling, 
the principle limitation of this approach is the exponen- 
tial decay of (SWAP^) with boundary length. Following 
Ref. HU a "ratio trick" may be employed in which the 
ratio of expectation values of SWAP can be computed 
more efficiently: 



(SWAP^) = , Sa _ Sa 
(SWAP a ) 



x(L A ,-L A ) 



(11) 



where region A' is larger than A but with perimeter La> 
sufficiently close to La such that the ratio is not too 
small. While a single ratio combination gives the dif- 
ference of the entanglement entropy of two regions, the 
entanglement entropy of a single (large) region, required 
to compute the TEE, can be determined by computing 
the entropy of a single small region and adding succes- 
sively the computation of differences (swap ratios): (i.e. 
S An = Sao + ELi^ - for A n > • • • > A, > 

•••> Ao). 

The ratio method can be specialized to the dimer sys- 
tem at generalized RK points. An example of two over- 
lapping regions A' and A is shown in Fig. [I] Because the 
estimator used to compute (SWAP^) (Eq. [lQ[) does not 
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have support everywhere in the doubled Hilbert space 
(it is non-zero only for swappable configurations), a diffi- 
culty with directly computing the ratio in Eq. [TT] is that 
swappable configurations over A' are not a subset of those 
over A or vice-versa. Instead, what is actually directly 
possible to compute are the following ratios: 

_ (SWAP A /SWAP A ) _ (SWAP A /SWAP A ) 

11 — 77777^ a v i R — 



(SWAP^ 



(SWAP A ) 



(12) 

From these, the ratio we want is simply R/R . Each of 
these ratios has a simple Monte Carlo interpretation. For 
example, inserting 1 = ^ c \C)(C\ into the numerator 
of R, and acting each SWAP towards the center readily 
leads to: 



R 



ZmA)S(C\ Af )e(- E ^- E ^)) 

c_ 

^ e -E(c A )-E(C) S (c\ A ) 
c 

= J2e- AE ^S(C\ AI )U A (C) 

C\ A 



(13) 



where the sums have been explicitly written as restricted 
to A-swappable configurations, AE(C\ A >) = E{C A i) — 
E(C), and 



U A (C) 



exp(-E(C A ) - E{C)) 

£ exp(-E(C A ) - E(C)Y 
C\ A 



(14) 



To compute R, we sample dimer configurations (weighted 
by e - E ( c A)-E(C)^ ^ n suc ] 1 a wa y that the A-boundary al- 
ways remains swappable, and then measure the estimator 
£(CV)e~ AjB(CU/) . At the standard RK point one may 
set exp — » 1, and R is simply the ratio of configura- 
tions that are swappable for both A and A! to those only 
A-swappable. To compute R \ the region which is kept 
swappable is reversed. 

To generate configurations which always remain swap- 
pable over A, updates are done independently for each 
copy, except for certain constrained plaquettes which 
must be flipped simultaneously in copy 1 and 2. The 
constrained plaquettes are those which are not entirely in 
region A or B (the complement of A) . Some examples are 
shown in Fig. [I] For simply-connected bipartitions these 
updates should be ergodic over all possible A-swappable 
configurations. However, if A is not simply connected, 
it can be shown that this is not be the case, and an- 
other method for ratio updates are needed. Details are 
presented in the supplementary material.^ 



III. COMPUTATION OF TOPOLOGICAL 
ENTANGLEMENT ENTROPY 

The TEE has been computed at the RK point in 
Ref. [12| and [13] using Kasteleyn matrices P^l Here we re- 
port computations of this quantity using Monte Carlo, 
which also allows us to sample generalized RK points 





FIG. 2. Sample parallelogram extrapolation (left); strips 
with periodic boundary conditions assumed (right). 





FIG. 3. Levin- Wen (left), Kitaev-Preskill (right) construc- 
tions. 



whose ground states exhibit a quantum phase transition. 
Our findings are useful for further numerical studies away 
from RK-like wavefunctions. 

There are two routes to compute the TEE term nu- 
merically. The first is to use extrapolations: that is, to 
extrapolate the linear part of the entanglement entropy 
scaling for different perimeter-sized bipartitions and de- 
duce the intercept. The second is to make use of can- 
cellations: i.e., to consider a difference of bipartitions 
whose total net perimeter cancels while the net number 
of boundaries do not, thereby leaving the topological con- 
tribution 7. 

Using extrapolations, one has the option of employing 
simply connected polygonal bipartitions such as parallel- 
ograms of various sizes, or splitting the torus-the topol- 
ogy of periodic boundary conditions-by a non-simply 
connected strip of different widths (Fig. [2|. The former 
turns out to suffer from corner contributions masking the 
TEE which we describe in detail in Sec. IIVI We found 
extrapolations from strips, which have smooth bound- 
aries, suffer from even-odd finite-size effects for thin 
tori, thereby frustrating accurate intercept fits. There- 
fore, cancellation strategies turn out to be more use- 
ful for Monte Carlo simulations. Two cancellation ge- 
ometries are typically used: a Levin- Wen type con- 
struction^ shown in Fig. [3] gives the TEE as two pairs 
-27 = (Sabcd ~ Sabc) + {Sac ~ S A dc), while the 
Kitaev-Preskill construction 2 (Fig. [3| gives the TEE as 
three pairs plus an extra region: —7 = (Sabc — Sab) + 
(Sa — Sac) + (Sb — Sbc) + Sc- It turns out that in both 
of these constructions the number and type of corners 
cancel as well. 

In practice, however, using the types of updates de- 
scribed in Sec. [ill the Levin- Wen construction turns out 
to not be possible for the ratio method. The issue is 
that the updates fail to be ergodic in the configuration 
space. The problem occurs for calculations when the re- 
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FIG. 4. (Color Online) Triangular lattice version of the 
Kitaev-Preskill construction used in the Monte Carlo runs. 



TABLE I. Numerical values for selected points of Figure [5] 





Monte Carlo Run | 


7/ In 2 


w 


= 6 (18 x 18 lattice) 


1.001 ±0.003 


w 


= 6 (16 x 16 lattice) 


0.994 ± 0.002 


w 


= 5 (18 x 18 lattice) 


0.995 ± 0.003 


w 


= 5 (16 x 16 lattice) 


0.995 ±0.001 
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FIG. 5. (Color Online) Net entanglement entropy (TEE) from 
the Kitaev-Preskill construction, for the RK wavefunction on 
the triangular lattice, computed using the partitions shown 
in figure [4] Unless noted, that ground state parity Q is (0, 0). 
The solid horizontal line is the predicted value of the TEE 
= In 2. 



gions Sabcd and Sac are the constrained regions (de- 
nominators in the ratios of Sec. [TTJ see supplementary 
material 2 ^). Consequently the Kitaev-Preskill construc- 
tion offers the best method for extraction of 7. 



A. RK point 

We take region ABC to be an inscribed regular 
hexagon with divisions as shown in Fig. [4j w is the outer 
hexagon side width. Our results for the RK point wave- 
function are summarized in Fig. [5] and Table [T] for vari- 
ous sized hexagons and lattices. If the spacing from the 
outer hexagon edge to the lattice edge is kept on order 
of the hexagon side width or greater, the resulting TEE 
depends little on lattice sizes, except in limiting the size 
of the outer hexagon that can be used (this turns out 
not to be the case for more general wavefunctions) . For 
w = 6 convergence to the expected value of In 2 is seen 
(within the error bars). We also computed the TEE for 
the different ground sta tes o n the torus, ^ = (1,0), (1,1) 
(in the notation of Sec. I A). These states also converge 
to the expected value, In 2. 



B. TEE approaching a quantum phase transition: 
generalized RK wavefunction 

A number of p hases fr om other classical weightings 
have been studied ! 16 * 19 * 20 ^ One class of modifications in- 
terpolate between the liquid-like RK point and sym- 
metry broken phases by preferentially weighting dimers 
on links reflecting the desired ordered configuration. 20 
These wavefunctions, however, explicitly break transla- 
tion invariance and therefore constructions which rely on 



perimeter cancellations fail to accurately probe the TEE 
unless very large regions can be used. 

A classical dimer model that preserves translation in- 
variance was discussed by Trousselet et al. In this model 
the wavefunction was weighted to favor dimers in paral- 
lel (equivalent to a flippable plaquette), as if they would 
"interact classically" in the sense of E(C) as a classical 
energy. In this case, E(C) = ln(a)iV/(C), where Nf(C) 
is the number of flippable plaquettes for a covering C, 
and a is an adjustable parameter, a = 1 corresponds 
to the RK point while a < 1 favors flippable plaquettes. 
Such interacting wav efunct ions are the groundstate of the 
RK-like HamiltonianPESl 



H a = t y %2-(\&)(&\ + h.c.) 



(a-^// 2 |£7)(ZS7| +a AN f/ 2 \/S?)(/S/\) 



(15) 



with ANf = Nf(&)-Nf(&) (the difference in flippable 
plaquettes due to a flip of plaquette p). 

The topological liquid phase with a finite correlation 
length persists for a < 1 until the wavefunction under- 
goes a first order phase transition near a c ~ 0.2! 16 * 2Q I 
Below this value only configurations that have the max- 
imum number of flippable plaquettes contribute to the 
wavefunction. These are 12 columnar symmetry-broken 
configurations, plus a large number of configurations re- 
lated by single line shifts across the lattice EHUD 

Using the ratio method and the Kitaev-Preskill con- 
struction we are able to probe the TEE of such inter- 
acting wavefunctions in the liquid regime as the system 
approaches the transition point. Although the TEE is 
believed to be universal in the topological liquid regime, 
however, the only known example at T = was given 
by Stephan et alW These authors computed the TEE for 
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FIG. 6. (Color Online) TEE for the ground state of Eq.[l5|as 
a function of the parameter a. All points are computed on an 
18 x 18 lattice, with the TEE extracted from a Kitaev-Preskill 
construction with outer hexagon width w. 
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FIG. 7. (Color Online) Negative log of the dimer-dimer corre- 
lation at 6 lattice spacings apart on an 18 x 18 lattice for the 
ground state of Eq. 15 as a function of the parameter a. ||, _L 



refer to dimers relatively parallel or offset by 60° respectively. 



a dimer model starting on the triangular lattice (at the 
RK point) and interpolated to the square (critical) lat- 
tice using Kasteleyn matrices ! 29 * 30 ^ At a critical point the 
entanglement entropy is predicted to have a constant pos- 
itive shift. However Stephan et al. found that before ris- 
ing from the value of — In 2 for the triangular lattice at 
the RK point, the TEE first deviated from this value by 
decreasing. The flow was well described by a single com- 
bination of parameters t ■ L, where L is simultaneously 
the cylinder circumference (the geometry used) and bi- 
partition boundary, while t is the fugacity for dimers on 
diagonal bonds (t = is the square lattice fugacity; note 
this is a different t then in Eq. [T]) . Deviation away from 
In 2 begins for as large a value as t • L ~ 9 (for Renyi 
parameter n = 1.5). 

Our results for the interacting wavefunction are pre- 
sented in Fig. [6] as a function of the parameter a. Since 
the ordered state involves many local minima, we can- 
not go through the transition without non-local updates, 
which we have not implemented in our variant of the ra- 
tio method. The TEE appears to be a robust indicator of 
topological order for a reasonable range of a. However, 
below a ~ .5 stronger finite size effects begin to affect 
the result. An investigation of dimer-dimer correlations, 
discussed below, suggests that larger correlations are the 
main reason for the deviation. While the scaling with 
outer hexagon width w is not clear for the system sizes 
studied (inset to Fig. [6|, it suggests convergence to the 
expected value of In 2. 

The value where the TEE begins to diverge from the 
theoretical value, occurs when the bipartite region length 
scale, w, is on order 10 x the correlation length £ (for the 
current Renyi index n = 2). For the interacting wave- 
function (ground state of Eq. [l5| , the dimer-dimer cor- 
relation length does not diverge, since the transition is 
first order, however, it grows compared to the RK val- 
ues. Figure [7] shows the negative log of dimer-dimer cor- 



relation at 6 lattice spaces apart-the length scale of the 
boundary bipartitions w on an 18 x 18 lattice (the lattice 
size for the results on Fig. [6|. For a lattice of this size, 
dimer-dimer correlations do not converge to a perfect ex- 
ponential. Therefore the negative log of the correlation 
value at a characteristic length, w(= 6), is used to esti- 
mate w/£, where £ is an effective correlation length. The 
divergence from — In 2, beginning around a = .5, corre- 
sponds to a value of 8 in Fig. [7[ Expressed in this way, 
the value can be compared to the results of Stephan et 
al p where a deviation can also be seen at similar values 
(since t ~ l/£, 31 and deviations begin near L-t ~ O(10)). 
Taken together, these results suggest that at least for 
dimer systems, the TEE is quite sensitive to finite cor- 
relations, requiring bipartition length scales up to O(10) 
times the correlation length. If such scaling holds more 
generally, it would constitute an important limitation on 
numerics to extract the TEE. 



IV. CORNER CONTRIBUTIONS TO 
ENTANGLEMENT ENTROPY 

As noted earlier, bipartitions with corners contribute a 
non-universal constant, ft, to the entanglement entropy. 
The presence of corner shifts can be readily seen in the 
scaling of various polygonal regions, in particular, of the 
hexagons, triangles, and rhombi shown in Fig. [8j The 
offset appears to be a constant shift. This is more pre- 
cisely measured by computing differences of bipartitions 
as discussed below. 

Corner effects have been studied for the integer quan- 
tum hall wavefunction Pin other topological phases there 
have not been any studies on corner contributions that 
we are aware of, although Ref. [13| also noted the poten- 
tial for non-universal constant corner contributions away 
from the critical point. Our results are consistent with a 
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FIG. 9. (Color Online) Six lattice corners. Type I corner 
(top), left to right: a, f3 a. Corresponding type II on bottom. 



FIG. 8. (Color Online) Entanglement entropy and best fit 
line as a function of polygonal region perimeters, L, show- 
ing offsets for different polygonal regions. Fits include total 
constant term; for comparison the expected 7 = In 2 = .693 
(solid line in the inset). 



TABLE II. Corner shifts a%. The estimate is computed from 
the analysis described in the appendix. See text and figure [To] 
for details on Monte Carlo runs employed. 



Cornerl Estimate 



Measured 



la 
lla 
1/3 
110 
la 
lla 



-0.02 
0.09 
-0.26 
-0.02 
-0.62 
-0.23 



0.008 
0.0946 
-0.2592 
-0.0011 
-0.677 
-0.238 



0.002 

0.0001 

0.0004 

0.0004 

0.002 

0.002 



total contribution k of the form k = 

QjiTli as found in 

i 

Ref. [9j where i labels the types of corner, rii the number 
of corners of type i while {a^} are coefficients to be deter- 
mined. With this form, k exactly cancels for the Kitaev- 
Preskill construction if, as required by the symmetry of 
entanglement entropy ai = where i is the comple- 
ment of i. Apart from the numerical values, the most 
important result of this section, is that corner effects do 
saturate to constant shifts, allowing for the extraction of 
the TEE by the approach of canceling regions. 

Since the linear shifts in Fig. [5] include the topologi- 
cal term, the triangles actually have the smallest shift, 
despite having sharper corners. More precise results con- 
firm this. To see why sharper corners can have a small 
effect, first note that unlike the continuum, the types 
of corners i are not solely determined by angle. Micro- 
scopically, on the triangular lattice there are a total of 
6 types of corners (plus complements) labeled as i G la, 
1/3, la, lla, II/3, and lla. The Roman numeral refers 
to the angle I = 60°, II = 120° and the Greek let- 
ter as to whether incoming links are included or not, 
as shown in Fig. [9j Measured values of the shifts ai 
are presented in Table [ill Estimates are obtained from 



a mean-field-like analysis presented in the appendix and 
supplementary material!^ Two important points can be 
deduced from the analysis that give intuition both as to 
the form of k and for the original question, i.e. why 
sharp corners have a small contribution. First, the cor- 
ner terms can be thought of most naturally as part of 
a generalized linear scaling with boundary vertices {^}, 
Sa ~ J2j — J2j for j G {v} where i is the lat- 
tice spacing. When all ctj are equivalent (ctj = a), as 
in a straight side, this becomes a linear scaling = Laol. 
Then the corner shift can be understood as single dislo- 
cation a — » oti, so that ai = Aa = oti — a. Secondly, 
to first order, the scaling ctj is determined by the num- 
ber of links in region A versus B (the complement) that 
meet at boundary vertices. For side vertices away from 
corners these are 4 and 2, therefore corners la and II/3 
which have the same incoming links have the smallest 
shift, followed by (3, 3) which adds entropy for lla, then 
(1,5) for 1/3 and ILr, and (0,6) for la. These values are 
all in agreement with the ordering shown in Table [ill 

Several Monte Carlo runs were used to extract the mea- 
sured di. The types I (II) a were taken from combining 
the information from: 1) the linear scalings in Fig. [8] 
removing the TEE entropy (assuming it is In 2), 2) the 
differences between those polygons (note any combina- 
tion of those only gives 2 a (lla) — a(Ia)), and 3) specially 
constructed combinations called "la, lla combinations" 
shown in Fig. 
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The measured are easily ob- 

tained solely from the L-shaped combination shown in 
Fig. 10 ("1/3, II/3 combinations"); similarly for I(II)a. 

From Fig. [To) the most important feature is that corner 
effects do indeed saturate to constant shifts, rendering 
constructions like the Kitaev-Preskill useful for extract 
the TEE. 



V. CONCLUSIONS 

We have demonstrated that Monte Carlo techniques 
provide a viable method for the numerical calculation of 
the TEE of the ground states of generalized RK points. 
In addition to confirming existing results, we have also 
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FIG. 10. (Color Online) Bipartition geometries used to extract corner contributions. Entropy differences are computed between 
the black and red regions. The dotted line represent sides with "missing links" used to form f3 corners. The la and Ila 
combinations involve pairs plus an extra region, thereby requiring a subtraction of the TEE (taken to be In 2). The blue 
horizontal lines are the values used for table UTI 



been able to apply the method to a generalized RK wave- 
function and thereby investigate the behavior of TEE ap- 
proaching a quantum phase transition. Our results sug- 
gest that the TEE is indeed a robust indicator of topo- 
logical order in the thermodynamic limit throughout the 
liquid regime. However, we also find a strong depen- 
dence on correlations that requires bipartitions with side 
lengths of order 10 times the correlation length. If this 
estimate applies to other systems, this implies that nu- 
merics will be severely constrained in the vicinity of a 
second order phase transition. These results are there- 
fore an important guide for future quantum Monte Carlo 
studies. 

The third aspect that we have been able to study here 
is the potential for corner contributions. First the magni- 
tude of corner shifts are of order 7 so their effects must be 
controlled either by canceling them out or by dealing ex- 
clusively with smooth boundaries. On a lattice with pe- 
riodic boundary conditions, however, these are not topo- 
logically trivial. For example, difficulties in the linear ex- 
trapolations of Ref.[T2lmav have been due to the effective 
variation of edges associated with the radial-like region 
definitions. Second, at least for the RK point, we can 
verify that corner shifts do indeed saturate to constant 
values, thereby making possible constructions which can- 
cel corners linearly. It remains to be seen whether this is 
the case for more general wavefunctions. 

Our results are important for future work employing 
quantum Monte Carlo and for conclusions about the 
practical utility for TEE as probe of topological order. 



In this regard, the present results allow us to conclude 
that TEE is a useful tool as long as correlation lengths 
are small enough in at least one point of the topologi- 
cally ordered liquid phase. While the only other known 
example of the behavior of the TEE on approaching a 
quantum phase transition does suggest a similar behav- 
ior,^ further numerical studies of TEE in exotic phases 
of matter should be of great value. 
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Appendix A: Estimate of corner effects. 

Using Kasteleyn matrices] 29 * 30 * an estimate for corner 
effects can obtained. Kasteleyn showed that the number 
of dimer coverings on a lattice is given by the Pfaffian of a 
matrix Kij, with z, j labeling lattice sites and ±1 (or 
suitable weights in generalizations) for nearest neighbors 
and zero otherwise; the sign is determined according to 
a convention of arrows placed on the lattice. We refer 
readers to Refs. [13] and [3TJ Only two facts will be used 
for the following arguments. First, Z = the number of 
coverings = the Pfaffian of K (written Pf(iT)). Second, 
the square of the Pfaffian is the determinant for an even- 
dimensional matrix: (Pf(i\")) 2 = Det(K). 




FIG. 11. (Color Online) Boundary vertices in black shown 
taking values {v} = {baabba} (periodic boundary conditions 
are not assumed). The lattice to be included in K — E are 
the darkened links. 



For a configuration to be swappable, one simply needs 
to satisfy the dimer constraint along boundary vertices 
= {v}, which are vertices that have links belonging to 
both B and A. For each boundary vertex in each config- 
uration, a dimer can be on a B link or an A link, which we 
will consider as a "side parity" labeled as either v = a,b 
if the vertex has a dimer in A or B. Then swappabil- 
ity simply imposes that for every boundary vertex, the 
dimer in copy 1 is in the same side as the dimer in copy 
2 (v\ = V2). Using a similar notation to Ref. [13j the 
number of swappable configurations can be written as: 



^SWAP 



£ 

{v} 



Z\\{ v } • Z 2 \ 



(Al) 



where Zi^\ v is the set of dimer coverings in copy 1(2) 
given a set of A, B occupations, or side parities of {v}; 
the sum runs over the combinations specifying a or b 
for the number of boundary vertices = N v . Note that 



Zi\ 



{v} 



J 2\{v} 



and each can be written as a Pfaffian, 



with certain constrained links removed. The links to be 
removed are those which touch a vertex v on the side 
A if v = 6, or on the side B if v = a (the Pfaffian will 
then generate combinations on the occupied side). An 
example is shown in Fig. 11 An "exclusion matrix" E 
can be defined so that for a removed link from site r to 
8, E rs = K rs and it is zero otherwise, so that: 

^swap = £(Pf(# " ^Im)) 2 = E Det (^ " ^Im) 
W M 

(A2) 

with K and E now referring to only a single copy. Finally 



from Eq. 10 the entanglement entropy can be written as: 
■Zswap 

= — m 

{«} 



S A = -ln- 



= -*L, Det{K) (A3) 



which matches the n = 2 Renyi entropy in Ref. 13. Next 
we can employ a perturbation theory of matrice^' and 



the trace identity, 

Bet(K-E\ {v} ) _ DetOFQDettl-if- 1 ^}) 
Det(K) ~~ Det(if) 
= exp(Trln(l-if- 1 i?| w )) 

= exp(Tr(- J FC"- 1 £: - -K^EK^E ...)), 



(A4) 



with E\{ v y implied in the last line. K 1 can be diagonal- 
ized in Fourier space It decays exponentially in vertex 
separation; for nearest neighbors z, j K~j = ±1/6, it van- 
ishes for next to nearest, and is ~ 0.02 for a three-link 
separation. Because of the number of terms contributing 
to the trace also grows, the series does not converge as 
fast as K~ x ~ 1/6, but roughly ~ 1/2, so that at least 
these two terms should be kept, and we expect the es- 
timates to be accurate to the order of 10% if we ignore 
higher order terms of the expansion. Ignoring those terms 
and making the definitions Ti|^} =Tr(— K~ x E\{ v }) and 
T 2 | w = Tr(-lK- 1 E\ {v} K- 1 E\ {v} ), the resulting form 
is 



5A = -ln^exp (Ti| 



{v} 



T 2\{v})- 



(A5) 



{v} 



The sum ranges over the a, b choices for each vertex. 

If Ti(2) depended only on the number of a versus b 
vertices Eq. |A5| could be written as a binomial expan- 
sion or Gaussian in the large N v limit. However, in gen- 
eral Ti(2) depend on the details (ordering) of a particu- 
lar set of a, b choices of {v}. Specifically, the first term 
T\ = (K~ 1 )ij(E\{ v y)ji is only non-zero for z, j nearest 
neighbors of the "excluded lattice". It is evaluated as 
~\^li z i — ^1? where the sum is over all vertices z, 
Z{ is the number of nearest neighbors for the zth ver- 
tex (of the excluded lattice). T\ can readily summed to 
— | x (number of removed links). This value depends on 
the number of a versus b values and also the ordering. 
For example, the number of links in the "excluded lat- 
tice" increases by 2 for every v = 6, but may increase by 
3 or 4 for v = a, depending on whether it follows an a 
or 6, respectively. Therefore to proceed, one can define 
a mean-field contribution as — | x 2 per v = b vertex 
while — | x (3P + 4(1 — P)) per v = a vertex where P 
is the likelihood of finding au = a neighbor, which can 
be determined self-consistently. The second term T 2 can 
be put in a similar form provided one ignores K~ x terms 
connecting three links or more, as this contains a factor of 
0.02. Then T 2 can be shown to be — ^ XlJ^f + — !)]• 
One can also make a mean field estimate for the average 
contribution to T 2 by considering combinations of neigh- 
boring preceding vertices. Note, that the form of these 
terms as a sum suggests the pot enti al for the generalized 
linear scaling mentioned in Sec. |IV[ Sa ~ o^di- 

If we call the average contribution to T per v = a 
vertex on a straight side, = r s id e ,a and similarly per v = b 
vertex = r s ^e,6? then given a combination of vertices {v} 
(along a side) having n a , v = a vertices, T = T\ +T 2 can 
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be written in terms of these average values: a's, n a , in {v}), the sum (Eq. |A5[ ) can be written as a 

binomial expansion and approximated as the Gaussian 
T = n a r S i de}a + (N v - n a )r side}b . (A6) integral in the large N v limit: 

Put into this form (depending only on the number of 
I 



r l~^2~ ( 2(n a - N v /2) 2 \ 
S A ~ - In / dn a J — — exp ( — h iV v ln(2) + r side ^(N v - n a ) + r side , a n a j 

(A7) 



Ar / i /ri\ 7~side,b H~ 'T side. a (j~side,a ^~side,b) \ / at- \ 

=iv v - ln(2) cta(A/ v ), 



Finally a corner shift for a corner of type c, can be ex- 
tracted from the difference ASa(c) = Sa(N v — 1 side ver- 
tices +c ) —Sa(N v side vertices); with the exception of la 
discussed below. One can think of the first term as a bi- 
nomial (Gaussian) sum over N v — 1 non-corner terms mul- 
tiplied by a factor of x(e Tc a + e Tc ' b ) which generates the 
last two possibilities at the corner with similar mean-field 
contributions for corner vertices: T\ a ^ ai tiq,^, ri^ ?a , r^^, 
etc!^ Then we can write: 

AS A = <t(N v - 1) - ln(e Tc - + e T ^) - a(N v ), (A8) 



or, 



AS A = ln(2) 



7~side,b H~ ^~side,a {^~side,a ^~side,b) 



-ln(e r ^ +e Tc - b ). 



(A9) 



For the case of la, since at the "corner" there is no bound- 



ary vertex. To compute the shift in this case, instead two 
vertices are removed, and the factor r^ a a refers to a sub- 
stitution of the vertex just before the corner. So the 
above formulae can be used with r*'s and the first term 
(enclosed in parentheses) multiplied by 2. 



Inserting estimates for the r's (see supplementary ma- 
terial) yields the results given in Tab. [n] and are m 
good agreement with the measured valued. Also, note 
that the term multiplying N v in cta(N v ) (Eq. |A7| ) is an 
estimation of the slope of the entanglement entropy. The 
value for the slope obtained by plugging in the mean-field 
values gives 0.51, which is on order 10% from the fitted 
slope values near 0.58. The relative success of the sub- 
stitutions suggest that the corners can be seen explicitly 
in the generalized linear scaling J2i a i^i noted above and 
in Sec. HV1 
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I. MEAN-FIELD APPROXIMATIONS FOR 
CORNER EFFECTS 

We can estimate on average the contribution to the 
traces Ti, T 2 , per vertex type r. The results in terms 
of the probability for finding a vertex with side parity 
v = a, P, are shown in Tab. I. Each row in the ta- 
ble represents the average incremental contribution to T 
(= Ti -f T 2 ) for a particular vertex, say x-a label for every 
row which includes the type of vertex, and whether it is 
occupied a or b (i.e. x = "side a, 6"). Then for each x, 
the coefficients for the terms involving the 1/6 factor are 
the determined as increment in the value of — 
(i.e. Ti) by computing the increase to T\ adding an x. 
Specifically, we start with a "single vertex exclusion lat- 
tice" (a single inner vertex, its links, and neighboring 
vertices) which can be comprised of an a or b vertex (oc- 
curring with probabilities P and (1 — P)). For each case 
the increment in T\ is computed after adding a second 
vertex x to make xa and xb respectively. Each term 
enters with probability P and 1 — P, thus giving the 
first two coefficients in the table (or a single number if 
the coefficients are the same). Similarly the coefficients 
for the second order terms are the increase of the sum 
(-1/72) Y,i zf + z%(zi - 1) (i-e. T 2 ) by comparing a two- 
vertex exclusion lattice comprised of aa, aft, 6a, bb (occur- 
ring with probabilities P 2 , P(l - P), P(l - P), (1 - P) 2 ), 
and adding a third vertex x to make xaa, xab, xba, and 
xbb respectively. When x is a corner that is asymmetric, 
such as f3 corners, an average (equal weight) over each 
side is also taken for each term. 

Finally to obtain, P, we note that before the inte- 
gration in Eq. A7, the Gaussian has mean \i = PN V1 
which defines the mean-field probability P used to self- 
consistently determine the r's to be: 

p 1 _j_ ^~side,a{P^) ^~side,b{P^) 

2 4 

With the estimates of Tab. I, numerically, P is: 

P = J (l29 - a/16369) « 0.2647. (2) 

II. CONFIGURATIONS WHICH CANNOT BE 
REACHED BY THE RATIO METHOD UPDATES 

To see that there are configurations which cannot be 
reached (for non-simply-connected regions) with our im- 



TABLE I. Mean- field estimate of the average contribution to 
the trace T per vertex of type t with a value a or b respectively 
= Tt,a{b)' P is the likelihood of finding an a- valued vertex 
determined self-consistently (Eqs. 1 and 2). r^ a a i s computed 
for the vertex just before the corner. See text. 



plementation of the ratio method, it is easiest to think 
of the transition loop picture between configurations. 1 3 
Starting from a reference configuration, any other config- 
uration may be described by flipping dimer occupations 
along closed loops (transition loops). Defining an oper- 
ator (Tij which flips the dimer occupation on a link con- 
necting i to j, transition loops i = {i,j} where {i, j} are 
a set of links forming closed loops, can be found connect- 
ing any two configurations: 

ic>=n^ic>- (3) 

£ 

In this way, dimer coverings are seen as superpositions of 
loops over a reference configuration. Note that a plaque- 
tte flip can be considered a minimal (four-link) transition 
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loop. 



Now consider a region A comprising two-disconnected 
parts, which constitutes a denominator or constrained 
region of the ratio method. Figure 1 shows illustrations 
of a doubled-dimer system ket, with such a disconnected 
region, and various plaquette flips as transition loops in 
blue. For such a bipartition, one can define a line in 
each copy connecting the two regions shown in the fig- 
ure by the dashed red lines. The Monte Carlo updates 
that affect the number of transition-loop-crossings over 
these lines, is such that the sum of crossings over both 
copies can only change by two. To show this, Fig. 1 illus- 
trates the three unique types of moves that could affect 
the crossings. In particular, because of the constrained 
plaquettes, the bottom panel includes plaquette flips that 
must be simultaneously applied to each copy. 



FIG. 1. (Color Online) A tensor product of copy 1 (left) 
and 2 (right), with transition loops for three types of Monte 
Carlo updates (top, middle, bottom) which affect the parity 
of transition-loop-crossings through the red dashed line. 




However there do exist configurations which are both 
swappable but contain only an odd number of crossings. 
One example is shown in Fig. 2. Furthermore, these con- 
figurations could be reached if plaquette flips were un- 
constrained and therefore should be part of the ground 
state. Importantly there is no such conserved parity if 
the second disconnected piece of A were not present. In 
that case any loop surrounding a region can be removed 
by outward plaquette flips and contracted to the identity 
in each copy independently. 



FIG. 2. (Color Online) A tensor product of copy 1 (left) 
and 2 (right) with a transition loop in blue which cannot be 
generated by the ratio method but which should belong to 
the (doubled) ground state. See text for explanation. 



These sorts of regions occur in the Levin- Wen construc- 
tion, 4 and are therefore not amenable to our sampling via 
the ratio method without non-local updates. 
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